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Density functional theory (DFT) calculations of 58 liquid elements at their triple point show that 
most metals exhibit near proportionality between thermal fluctuations between virial and potential- 
energy in the isochoric ensemble. This demonstrates a general “hidden” scale invariance of metals 
making the dense part of the thermodynamic phase diagram effectively one dimensional with respect 
to structure and dynamics. DFT computed density scaling exponents, related to the Griineisen pa¬ 
rameter, are in good agreement with experimental values for 16 elements where reliable data were 
available. Hidden scale invariance is demonstrated in detail for magnesium by showing invariance 
of structure and dynamics. Computed melting curves of period three metals follow curves with 
invariance (isomorphs). The experimental structure factor of magnesium is predicted by assuming 
scale invariant inverse power-law (IPL) pair interactions. However, crystal packings of several tran¬ 
sition metals (V, Cr, Mn, Fe, Nb, Mo, Ta, W and Hg), most post-transition metals (Ga, In, Sn, 
and Tl) and the metalloids Si and Ge cannot be explained by the IPL assumption. Thus, hidden 
scale invariance can be present even when the IPL-approximation is inadequate. The virial-energy 
correlation coefficient of iron and phosphorous is shown to increase at elevated pressures. Finally, we 
discuss how scale invariance explains the Griineisen equation of state and a number of well-known 
empirical melting and freezing rules. 


I. INTRODUCTION 

Scale invariance plays an important role in many 
branches of science. It greatly simplifies a given phe¬ 
nomenon by reducing the parameter space and introduc¬ 
ing universalities over length or time scales. Some ex¬ 
amples of this are the size distribution of earth quakes 
[T], Brownian motions of microscopic particles [5], cos¬ 
mic microwave background radiation Em, and biological 
fractal structures E] such as those of lung tissue or Ro- 
manesco broccoli. In the physics of matter, scale invari¬ 
ance controls the properties of a fluid near the gas-liquid 
critical point El- This paper establishes an approximate 
“hidden” scale invariance in the dense liquid part of the 
thermodynamic phase diagram of certain elements from 
ab initio computations. While the above classic examples 
of scale invariance cover several decades of variations, the 
scale invariance of condensed matter covers much smaller 
length-scale differences. Nonetheless, hidden scale invari¬ 
ance implies that the phase diagram becomes effectively 
one-dimensional. 

Specifically, we have performed ab initio density func¬ 
tional theory (DFT) computations on 58 liquid elements 
at their triple point. We infer hidden scale invariance of 
metals from strong correlations in thermal fluctuations 
of virial W (the potential part of the pressure) and in¬ 
ternal potential energy U. We find that metals in gen¬ 
eral posses hidden scale invariance. These results give 
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the first ab initio quantum-mechanical confirmation of 
the picture proposed in Refs. EHia according to which 
some systems have strong virial potential-energy corre¬ 
lations in the condensed phases, whereas this is not the 
case for other systems - typically those with strong di¬ 
rectional bonding. Liquids belonging to the former class 
of systems were originally referred to as “strongly corre¬ 
lating”, but to avoid confusion with strongly correlated 
quantum systems the term Roskilde (R) simple system 
is now sometimes used [HHin]. Although the focus here 
is on monatomic systems, experimental and molecular- 
dynamics simulation results have shown that (R) simple 
systems include some molecular van der Waals bonded 
systems [5D], polymers [niiia, and crystalline solids [53] . 
Returning to the elements, molecular liquids like N 2 Ej 
and the noble gases [7| are also expected to be simple |S] . 

What is the possible microscopic origin of hidden scale 
invariance? To answer this, we start by noting that at 
high pressure the dominant interatomic forces are harshly 
repulsive |5^H55] . These forces can be modeled approx¬ 
imately by scale-invariant inverse-power-law (IPL) pair- 
potentials terms oc r”" [33M3] in which r is the inter¬ 
particle distance, plus a density-dependent constant g{p) 
taking long-ranged attractive interactions into account 
ElETj- By Euler’s theorem for homogeneous functions, 
there is exact correlation between the fluctuations of 
virial and potential-energy in the IPL case. Strong corre¬ 
lations between virial and potential energy are, however, 
not necessarily a consequence of approximately scale- 
invariant pair interactions, but rather a property of the 
multidimensional energy-function U (R, V) per se. This 
fact motivates the present paper in which we show from 
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first-principles computations that many metals obey Eq. 
([^ also at low pressure, i.e., when virtually uncompressed 
compared to the zero-pressure state. 

For metals, the outer electrons result in complex many- 
body interactions, which this study takes explicitly into 
account. At low pressure hidden scale invariance is non¬ 
trivial since the forces on the atoms are not at all domi¬ 
nated by pairwise additive repulsive IPL-type forces |49j . 
Hidden scale invariance have been found to be absent 
if the multidimensional energy function 17(R, E) have 
prominent contributions characterized by two or more 
length scales 0 0. Examples of this are water with 
hydrogen bonds and core repulsions, and the Dzugutov 
model with core repulsions supplemented by repulsive 
forces on second-nearest neighbors 0 0. However, it 
is possible to have multi-body interactions and hidden 
scale invariance. Below we do not make any assump¬ 
tions on the potential-energy landscape and verify hidden 
scale invariance for the majority of metals directly from 
a quantum-mechanical description. Density Functional 
Theory (DFT) currently provides the best means to do 
so [5D1153) . and although DFT can only approximate the 
true electronic behavior, it provides accurate predictions 
on a broad scale, such as crystal structures, densities and 
melting points. 

The remainder of the paper is organized as follows: 
Sec. |H] gives the theoretical background of hidden scale 
invariance and Sec. |HI| gives details on the ab initio DFT 
method. Results are presented in Sec. EYi and in Sec. 
[V]we relate our findings to Griineisen equation of state, 
and well-known empirical melting and freezing rules. 


II. THEORY: HIDDEN SCALE INVARIANCE 
OF CONDENSED MATTER 

Certain condensed-matter systems are characterized 
by a “hidden scale invariance” as reflected in the follow¬ 
ing approximate representation of the potential-energy 
function at a density p = N/V 

U{R)^h{pM-k)+Ng{p), ( 1 ) 

where the coordinates of the system’s N particles have 
been merged into a single vector, R = (ri, r 2 ,..., r^r), 
and the reduced coordinate is R = p^/^R. The inten¬ 
sive functions h{p) and g{p) both have dimension energy 
and $ is a dimensionless, state-point-independent func¬ 
tion of the dimensionless variable R, i.e., a function that 
involves no lengths or energies. Physically, Eq. Q means 
that a change of density to a good approximation leads 
to a linear affine transformation of the high-dimensional 
potential-energy surface. Thus if temperature is adjusted 
in proportion to h{p), state points in the thermodynamic 
phase diagram are arrived at for which the molecules ac¬ 
cording to Newton’s laws move in the same way, except 
for a uniform scaling of space and time. Sets of such 
state points are referred to as isomorphs and along the 


isomorphs structure and dynamics are identical in prop¬ 
erly reduced units to a good approximation [10) . The 
local slope d\n{T)/d\n{p) of the isomorph is given by 

7 (p) = d\n{h{p))/dln{p), (2) 

and is refered to as the density scaling exponent ffn- 
usi EHEa Thus, the phase diagram becomes 

effectively one-dimensional, and density and temperature 
merge into a single parameter. 

Hidden scale invariance is revealed in thermal fluctua¬ 
tions at a single state point: There are two contributions 
to the pressure, the ideal-gas pressure pid - a term that is 
always present and which only depends on the particles’ 
(atoms’) velocities - plus a term deriving from the inter¬ 
action between the particles. The latter is the so-called 
virial W, and it only depends on the particles’ coordi¬ 
nates. W is an extensive quantity of dimension energy 
and the general pressure relation is p = pid + WjV. Hid¬ 
den scale invariance dictates that fluctuations of virial 
and potential energy are strongly correlated in the NVT 
ensemble 00 E]: 

VF(R, V) ^ j{p)U(R, V) + const. (3) 

where the microscopic virial is defined by kF(R, F) = 
—VdU(R, V)jdV 1^1^. Below we use the Person cor¬ 
relation coefficient 

R = (AlT(R)AC/(R))/\/([AfF(R)]2)([A17(R)]2) (4) 

of virial and potential energy fluctuations to determine 
to what degree the multi body energy function C/(R) 
have hidden scale invariance. Here, (...) indicate a ther¬ 
modynamic average in the NVT ensemble and A’s are 
differences to average values. R is an intensive quantity 
with a value between -1 and 1. A value close to 1 indicate 
hidden scale invariance. The density scaling exponent is 
given by m 

7 = (A1F(R)AC7(R))/([AC7(R)]2). (5) 

We have previously conjectured that metals possess a 
hidden scale invariance, but this was based on assuming 
Lennard-Jones type pair potentials 0 . Ref. 0 showed 
via the Effective Medium Approximation that pure cop¬ 
per, as well as a magnesium alloy, exhibit hidden scale 
invariance. We do not make such assumptions in this 
paper, as detailed in the following section. 

HI. METHOD: DENSITY FUNCTIONAL 
THEORY 

DFT molecular-dynamics is a computationally efficient 
method with a quantum-mechanical treatment of the 
electron-density field. Examples of recent successes of 
DFT are computation of anharmonic contribution due 
to phonon-phonon interactions for fee crystals of metals 
m, dynamics of water dissociative chemisorption on a 
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Ni surface EH, and accurate computations of the melting 
points of period-three metals [62] . 

In this paper we present calculations of 58 elements us¬ 
ing the Vienna Ab-initio Simulation Package (VASP) [53] 
employing the projector augmented wave (PAW) method 
[64j with the frozen-core assumption and the Perdew- 
Burke-Enzerhof (PBE) exchange-correlation functional 
[55] , Dispersion corrections are not included, but these 
are not expected to be important for metals near the 
melting point |62] . We use periodic simulation cells con¬ 
taining 125 atoms for most elements except for the ele¬ 
ments Li, Na, Mg and A1 where we use 256, 250, 384 and 
256 atoms, respectively. Initial equilibration trajectories 
cover between 9 to 24 ps corresponding to several struc¬ 
tural relaxation times. Constant temperature is obtained 
with the Langevin thermostat with a coupling time of 
1 ps. Table m contains the electronic configurations of the 
calculated electrons of the employed potentials as well as 
the cutoff energy of the plane wave basis set. All calcu¬ 
lations were done non-spin polarized, and the Brillouin 
zone was sampled at its center the P point. Statistical 
uncertainties are estimated by dividing MD trajectories 
into statistically independent blocks. Details on estimat¬ 
ing DFT triple points are given in the Appendix. 


IV. RESULTS 


A. Correlated virial and potential energy 
fluctuations 


Fig.[T] shows the results of DFT computations on the 
first six period-three elements. Each subfigure gives a 
scatter plot of virial versus potential energy of configu¬ 
rations taken from NVT equilibrium simulations of the 
liquid phase at the triple points. For the metals (Na, 
Mg and Al) and the metalloid (Si) the scatter plots show 
strong correlations, implying that Eq. is obeyed to 
a good approximation. The value of the Pearson cor¬ 
relation coefficient R (Eq. (^) quantifies how well the 
approximation is obeyed mm- The non-metals P and 
S do not exhibit strong WU correlations. 


B. Invariant structure and dynamics 

To demonstrate hidden scale invariance directly, we 
investigated liquid magnesium in more detail. The re¬ 
sults are summarized in Fig. which studies isomor¬ 
phic as well as isochoric state points for temperatures 
between 900 K and 1600 K. Panels (a) and (b) give the 
radial distribution functions for the isomorphic and iso¬ 
choric state points, respectively, while (c) and (d) give the 
translational order parameter Q suggested by Truskett, 
Torquato and Debenedetti |7D| and the reduced diffusion 
constant. We see that structure and dynamics are al¬ 
most invariant along the isomorph, which confirms mag¬ 
nesium’s hidden scale invariance. 



Potential energy [eV] 

FIG. 1. Results from ab initio quantum-mechanical cal¬ 
culations of the first six period-three elements in the liquid 
phase at the triple point. Each subfigure shows a scatter plot 
of the potential energy U (R, V) and the corresponding virial 
1V(R, V) for configurations of the NVT ensemble. Fluctua¬ 
tions of virial and potential energy are strongly correlated for 
the metals Na, Mg and Al and the metalloid silicon Si, but 
not for the covalently bonded non-metals P and S. The strong 
correlation validates hidden scale invariance for metals. The 
subfigures report the Pearson correlation coefficient R (Eq. 
0 ) and the density-scaling exponent 7 of Eq. 0 determined 
as the linear-regression slope (Eq. 0). 





FIG. 2. Predicted isomorph invariants of magnesium along 
an isomorph and an isochore for the same temperature varia¬ 
tion. Panels (a) and (b) show results for the radial distribu¬ 
tion function g{r) in reduced units, panels (c) and (d) show 
the translational order parameter Q of Debenedetti et al. [201 
and the reduced diffusion constant D. The reduced unit struc¬ 
tural and dynamical quantities g{r), Q, and D all vary much 
less along the isomorph ((a) and (c)) than along the isochore 
((b) and (d)). 
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Element 

N 

Electrons 

Emax 

[eV] 

T [exp.] 
[K] 

P [exp.] 

[g/cm®] 

P 

[GPa] 

R 

7 

cr 

Ripl 

3 

Li 

256 

2s^ 

140.0 

470 [454] 

0.534 

[0.512] 

0.29 

0.69(0.08) 

1.0(0.1) 

1.92(0.28) 

0.71 

4 

Be 

125 

2s" 

308.8 

1850 [1560] 

1.570 

[1.690] 

-0.27 

0.75(0.03) 

1.8(0.3) 

1.16(0.09) 

0.50 

5 

B 

125 

2s"2p^ 

318.6 

2400 [2349] 

2.245 

[2.08] 

-1.69 

0.10(0.12) 

— 

— 

— 

6 

c 

125 

2s" 2p" 

273.9 

6000 [4800] 

1.514 

[1.37] 

-1.04 

-0.04(0.18) 

— 

— 

— 

11 

Na 

250 

3s^ 

62.1 

370 

[371] 

0.941 

[0.927] 

-0.02 

0.83(0.02) 

1.9(0.1) 

1.25(0.03) 

0.89 

12 

Mg 

384 

3s^ 

98.5 

900 

[923] 

1.581 

[1.584] 

0.07 

0.90(0.01) 

2.5(0.2) 

1.15(0.03) 

0.81 

13 

A1 

256 

3s" 3p^ 

116.4 

1000 

[933] 

2.352 

[2.375] 

-0.27 

0.86(0.06) 

4.0(0.4) 

1.02(0.02) 

0.90 

14 

Si 

125 

3s" 3p" 

245.3 

1700 

[1687] 

2.827 

[2.57] 

0.51 

0.79(0.03) 

4.1(0.2) 

0.94(0.02) 

0.79 

15 

P 

125 

3s" 3p" 

255.0 

650 

[317] 

1.850 

[1.74] 

-0.43 

0.01(0.17) 

— 

— 

— 

16 

S 

125 

3s"3p'‘ 

258.7 

550 

[388] 

1.784 

[1.819] 

0.07 

0.14(0.09) 

— 

— 

— 

19 

K 

125 

3p®4si 

116.7 

350 

[337] 

0.815 

[0.828] 

0.07 

0.84(0.12) 

1.6(0.3) 

1.45(0.27) 

0.84 

20 

Ca 

125 

3p®4s" 

266.6 

1200 

[1115] 

1.378 

[1.378] 

-0.44 

0.80(0.06) 

1.9(0.1) 

1.20(0.03) 

0.83 

21 

Sc 

125 

4s" 3d^ 

154.8 

1900 

[1814] 

2.800 

[2.80] 

-0.38 

0.63(0.24) 

1.4(0.5) 

1.43(0.56) 

0.81 

22 

Ti 

125 

4s" 3d" 

178.3 

3900 

[1941] 

4.161 

[4.11] 

-0.48 

0.78(0.08) 

2.0(0.2) 

0.96(0.04) 

0.88 

23 

V 

125 

3p®4s"d^ 

263.7 

2500 

[2183] 

5.773 

[5.5] 

-0.13 

0.81(0.07) 

2.6(0.6) 

1.02(0.12) 

0.81 

24 

Cr 

125 

3p'‘4sid'' 

265.7 

2600 

[2180] 

6.735 

[6.3] 

-1.55 

0.90(0.04) 

3.3(0.9) 

0.98(0.07) 

0.83 

25 

Mn 

125 

4s"3d5 

269.9 

2400 

[1519] 

7.973 

[5.95] 

0.66 

0.93(0.02) 

3.6(0.2) 

0.95(0.03) 

0.72 

26 

Fe 

125 

3p®4s"3d'’ 

267.9 

2400 

[1811] 

8.338 

[6.98] 

0.03 

0.95(0.02) 

3.6(0.1) 

1.00(0.01) 

0.90 

27 

Co 

125 

4s" 3d" 

268.0 

1870 

[1768] 

8.924 

[7.75] 

0.22 

0.93(0.01) 

3.5(0.1) 

1.06(0.01) 

0.94 

28 

Ni 

125 

4s" 3d® 

269.5 

2000 

[1728] 

8.189 

[7.81] 

-0.17 

0.92(0.03) 

3.5(0.3) 

1.03(0.02) 

0.96 

29 

Cu 

125 

4s"3di“ 

295.4 

1480 

[1358] 

8.020 

[8.02] 

-0.72 

0.90(0.02) 

4.1(0.2) 

1.01(0.01) 

0.94 

30 

Zn 

125 

4s"3d"“ 

276.7 

760 

[693] 

6.570 

[6.57] 

-0.61 

0.53(0.12) 

3.3(0.9) 

1.03(0.04) 

0.43 

31 

Ga 

125 

4s"4p" 

134.7 

500 

[303] 

5.967 

[6.095] 

-0.41 

0.74(0.09) 

3.3(0.5) 

1.09(0.04) 

0.65 

32 

Ge 

125 

4s"3d"“4p" 

310.3 

1250 

[1211] 

5.600 

[5.60] 

-1.46 

0.82(0.15) 

4.8(1.1) 

0.91(0.02) 

0.80 

33 

As 

125 

4s" 4p® 

208.7 

1300 

[1090] 

5.220 

[5.22] 

0.58 

-0.08(0.13) 

— 

— 

— 

34 

Se 

125 

4s^4p'^ 

211.6 

550 

[494] 

3.838 

[3.99] 

0.16 

-0.02(0.21) 

— 

— 

— 

37 

Rb 

125 

4p®5s" 

121.9 

340 

[312] 

1.460 

[1.46] 

0.01 

0.78(0.20) 

1.8(0.5) 

1.26(0.32) 

0.86 

38 

Sr 

125 

4s"4p®5s" 

229.4 

1100 

[1050] 

2.375 

[2.375] 

-0.24 

0.88(0.14) 

1.9(0.7) 

1.19(0.22) 

0.89 

39 

Y 

125 

4s"4p®5s"4d" 

202.6 

1850 

[1799] 

4.359 

[4.24] 

-0.47 

0.60(0.09) 

1.3(0.2) 

1.47(0.25) 

0.82 

40 

Zr 

125 

4s"4p®5s"4d" 

229.9 

2500 

[2128] 

6.305 

[5.8] 

0.01 

0.82(0.07) 

2.2(0.5) 

1.13(0.11) 

0.93 

41 

Nb 

125 

4p®5sffid'‘ 

208.6 

2900 

[2750] 

7.669 


0.96 

0.89(0.04) 

3.2(0.2) 

1.01(0.02) 

0.91 

42 

Mo 

125 

5s"4d® 

224.6 

3000 

[2896] 

9.330 

[9.33] 

0.42 

0.89(0.03) 

3.2(0.4) 

0.98(0.04) 

0.77 

43 

Tc 

125 

5s"4d® 

228.7 

2500 

[2430] 

10.606 


-0.43 

0.86(0.03) 

4.3(0.3) 

0.77(0.12) 

0.08 

44 

Ru 

125 

5s"4d" 

213.3 

2800 

[2607] 

11.200 

[10.65] 

1.18 

0.94(0.05) 

4.6(0.5) 

0.97(0.02) 

0.85 

45 

Rh 

125 

Ssffid® 

229.0 

2400 

[2237] 

10.807 

[10.7] 

-1.86 

0.88(0.07) 

5.2(0.4) 

0.96(0.01) 

0.89 

46 

Pd 

125 

Ssffid" 

250.9 

1900 

[1828] 

10.380 

[10.38] 

-1.32 

0.92(0.04) 

4.9(0.5) 

0.98(0.01) 

0.94 

47 

Ag 

125 

5sffid"“ 

249.8 

1350 

[1235] 

9.320 

[9.320] 

0.98 

0.89(0.04) 

4.7(0.4) 

1.00(0.01) 

0.96 

48 

Cd 

125 

5s"4di“ 

274.3 

650 

[594] 

7.996 

[7.996] 

0.61 

0.73(0.08) 

5.1(0.7) 

1.00(0.02) 

0.80 

49 

In 

125 

5s"5p" 

95.9 

600 

[430] 

6.859 

[7.02] 

0.28 

0.90(0.05) 

4.2(0.5) 

1.02(0.04) 

0.88 

50 

Sn 

125 

5s" 5p" 

103.2 

900 

[505] 

6.685 

[6.99] 

-0.89 

0.88(0.06) 

4.7(0.9) 

0.99(0.03) 

0.83 

51 

Sb 

125 

5s^5p^ 

172.1 

950 

[904] 

6.530 

[6.53] 

0.68 

0.40(0.15) 

3.0(1.7) 

1.03(0.33) 

0.41 

52 

Te 

125 

5s^5p‘* 

175.0 

750 

[723] 

5.700 

[5.70] 

-0.64 

0.14(0.30) 

— 

— 

— 

55 

Cs 

125 

5s^5p®6s^ 

220.3 

330 

[302] 

1.826 

[1.843] 

0.04 

0.90(0.22) 

1.7(0.6) 

1.33(0.77) 

0.96 

56 

Ba 

125 

5s^5p®6s^ 

187.2 

1050 

[1000] 

3.338 

[3.338] 

-0.41 

0.55(0.18) 

0.8(0.2) 

2.27(0.85) 

0.59 

57 

La 

125 

5s"5p®6s"5d" 

219.3 

1280 

[1193] 

5.940 

[5.94] 

-1.86 

0.72(0.23) 

1.7(0.6) 

1.25(0.30) 

0.86 

72 

Hf 

125 

5p'‘6s"5d" 

220.3 

2600 

[2506] 

12.349 

[12] 

-1.92 

0.83(0.03) 

2.2(0.3) 

1.16(0.07) 

0.91 

73 

Ta 

125 

5p'‘6s"5d® 

223.7 

3450 

[3290] 

15.000 

[15] 

0.45 

0.90(0.03) 

3.3(0.3) 

1.02(0.02) 

0.91 

74 

w 

125 

5p®6s"5d‘‘ 

223.1 

3900 

[3695] 

16.803 

[17.6] 

-0.78 

0.86(0.08) 

3.7(0.6) 

0.99(0.03) 

0.88 

75 

Re 

125 

6s" 5d® 

226.2 

3650 

[3459] 

18.787 

[18.9] 

0.98 

0.82(0.09) 

4.5(0.7) 

0.97(0.02) 

0.81 

76 

Os 

125 

6s" 5d® 

228.0 

3450 

[3306] 

19.741 

[20] 

0.46 

0.86(0.06) 

5.1(0.4) 

0.97(0.01) 

0.85 

77 

Ir 

125 

6s" 5d® 

210.9 

2900 

[2719] 

19.275 

[19] 

-0.91 

0.71(0.07) 

5.1(0.4) 

0.96(0.01) 

0.83 

78 

Pt 

125 

6s" 5d" 

230.3 

2200 

[2041] 

18.532 

[19.77] 

-1.00 

0.87(0.06) 

6.0(1.4) 

0.97(0.02) 

0.94 

79 

Au 

125 

6s"5d"“ 

229.9 

1470 

[1337] 

16.690 

[17.31] 

1.06 

0.86(0.14) 

7.9(1.6) 

0.95(0.02) 

0.92 

80 

Hg 

125 

6s"5d"“ 

233.2 

470 

[234] 

11.917 


0.48 

0.84(0.13) 

4.4(1.9) 

1.00(0.07) 

0.67 

81 

Tl 

125 

6s" 6p" 

90.1 

600 

[577] 

11.220 

[11.22] 

0.53 

0.90(0.09) 

3.7(0.8) 

1.06(0.04) 

0.86 

82 

Pb 

125 

6s" 6p" 

98.0 

900 

[601] 

10.671 

[10.66] 

0.48 

0.90(0.10) 

4.1(0.7) 

1.03(0.03) 

0.94 

83 

Bi 

125 

6s" 6p® 

105.0 

580 

[545] 

10.050 

[10.05] 

-0.27 

0.17(0.02) 

1.6(1.5) 

1.40(0.13) 

0.14 

84 

Po 

125 

6s" 6p^ 

159.7 

850 

[527] 

9.039 


-0.01 

0.36(0.21) 

2.5(1.7) 

1.13(0.40) 

0.35 


TABLE I. Correlation coefficients R (Eq. @) and density scaling exponents 7 (Eq. ([^) of the elements at the estimated DFT 
triple points T, p. Values in round parenthesis are the statistical uncertainties on R and 7 . N is the number of atoms in the 
periodic simulation cell. The electronic configuration of the calculated electrons as well as the plane wave cutoff energies Emax 
are listed. The last two columns gives the cr parameter of the IPL approximation (Eq. and the correlation between IPL and 
DFT energy fluctuations. Experimental values of T and p at the triple points are given in square brackets | 66H69| . 
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FIG. 3. Determination of the isomorphic temperatures at 
p/po = 1.3 for Mg using Eq. (j^ (with po = 1.581 g/cm^ 
being the triple point density). The ratio T/Tq is given as the 
ratio of the fluctuations of the inner energy of the scaled and 
unsealed configurations. 


tals generally have stronger virial potential-energy cor¬ 
relations than liquids [23]. We therefore conclude that 
metallic elements are (R) simple in the entire condensed- 
phase part of the phase diagram, i.e., exhibit hidden scale 
invariance. This excludes state points close to the criti¬ 
cal point, as well as those of the gas phase far from the 
melting line. 

Table |I| also reports the computed DFT scaling expo¬ 
nents 7 , which we discuss in details in the following Sec. 

IlYDl 

D. Density scaling exponents 

For a monatomic liquid of N classical particles (above 
the Debye temperature) the density scaling exponent 7 
can be determined purely from thermodynamic responses 

m- 


The isomorph in Fig is determined as follows: First 
we choose to identify state points with density increases 
of 9%, 20% and 30% relative to that of the triple point. 
Then temperatures along the isomorph are determined 
by relating the Boltzmann factors of scaled configurations 

m- 


exp{-U{Ro)/ksTo) oc exp(-C/(R)/fcBT). ( 6 ) 

For 100 representative configurations Rq at the reference 
point (i.e. the triple point) (To,po)j configurations are 
rescaled to the new density p and the energy of the scale 
configuration U{R) = C/(Ro[po/p]^^^) is evaluated with 
DFT. By taking the logarithm of Eq. (§ , we see that 
T/Tq is the ratio of the fluctuations of the energy at the 
reference point U (R) and the energy of the scaled config¬ 
urations C/(Ro). Fig.j^show an example of a scatter plot 
of these. The computed temperatures of the isomorphic 
states are 1110 K, 1370 K and 1600 K, respectively. 


C. The periodic table of hidden scale invariance 


The results of Fig. inspired us to study the elements 
in general to investigate whether all liquid metals ex¬ 
hibit hidden scale invariance. The results for 58 ele¬ 
ments are summarized in Fig. [^ and Table |l| We excluded 
some non-metallic elements (gray on Fig. for which 
standard semi-local density functionals are inaccurate. 
Metallic liquid elements all have strong or fairly strong 
virial potential-energy correlations at the triple point as 
quantified in the virial potential-energy correlation coef¬ 
ficient R. Most of the metals have a R larger than 80%, 
however, a few metallic elements show correlation in the 
range of 50%-60% (Li, Sc, Zn, Y, Ba). Scale invariance 
is expected to be worse for these elements. 

As mentioned, all systems have the hidden-scale- 
invariance property at high pressure where repulsive pair 
interactions dominate (see Sec. IVG). Moreover, crys- 


7 = [7G - kB/cy]/[l - 3kB/2cy] . (7) 

where 70 = ctpKx/pcy is the thermodynamic Griineisen 
parameter. Here, Up = [dV/dT]p/V is the isobaric ther¬ 
mal expansion coefficient, Kt = —V[dp/dV]T is the 
isothermal bulk modulus, c„ = Cy/N is the isochoric 
heat capacity per atom, and fcs is the Boltzmann con¬ 
stant. The density scaling exponent is proportional to 
the thermodynamic Griineisen parameter up to the ac¬ 
curacy of the Dulong-Petit approximation Cy ~ ZkB'- 
7 ~ — (see inset on Fig.[^. In Fig.j^we compare 

16 experimentally determined density scaling exponents 
m to the values computed with DFT for the liquid state 
(see Tab. |T]). Agreement is generally quite good. 

Incidentally, our computations show that the ac¬ 
claimed Lennard-Jones model does not properly reflect 
the physics of most metals because the DFT values of 
the density-scaling exponents are generally significantly 
smaller than 7 ~ 6 of the Lennard-Jones model |7]. 

Interestingly, with the exception of Fe and Tl, the el¬ 
ements Na, K, Rb and Cs forming the more open body 
centered cubic (bcc) crystal structure have lower scaling 
exponents than the elements Mg, Al, Co, Ni, Cu, Zn, Ag, 
Cd, Au and Pb that form close packed (cp) structures - 
face centered cubic (fee) or hexagonal close packed (hep). 
We will explore this further in the following Sec. |IVE[ 


E. Revisiting the inverse power-law model 

As mentioned in the introduction, hidden scale invari¬ 
ance may be explained if the interactions in a metal can 
be approximated with the IPL model [551 - H5] : U{R) ~ 
t/iPL(R) where 

N 

C/ipl(R) = Ng{p) + |r, - (8) 

i>j 

where A = kBTa^^p-^ is a materials dependent con¬ 
stant. The mean field term Ng{p) takes long-ranged 
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FIG. 4. Periodic table of elements with color coded baegrounds indicating the virial potential-energy correlation coefficient R. Greenish colors correspond to elements 
with strong correlations between fluctuations of virial and potential energy in the NVT ensemble of the liquid state at the triple point. These elements exhibit hidden 
scale invariance. The thick line separates metals and non-metals. Most metalloids and all non-metals investigated have weak to vanishing correlations. Gray elements 
were not considered since standard semi-local density functionals are inaccurate for these. See legend and Table [^for more details. 
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FIG. 5. Density scaling exponents from experimental ther¬ 
modynamic data collected in Ref. m compared to ab initio 
DFT computations for liquids. Elements forming bcc struc¬ 
tures are marked with blue diamonds, while elements forming 
one of the close-packed structures (fee or hep) are marked with 
green triangles. The inset compares the experimental 7 ’s to 
the Griineisen parameter 70 . The dashed line indicates the 
Dulong-Petit approximation {cv — 3fcs). 


attractive interactions into account El HZ]- Scale in¬ 
variance is not “hidden” for the IPL model, but a triv¬ 
ial result: Temperature and density merge into the di¬ 
mensionless parameter a = [Ap'^- i.e. a 
is the single parameter of the phase diagram and iso- 
morphs are given hy T (X p~^. To quantify the ac¬ 
curacy this approximation we list the correlation be¬ 
tween the IPL energy and the DFT energy, i?iPL = 
(AC/ipl(R)AC/(R))/V([AC/ipl(R)]2)([A17(R)]2), in Ta¬ 
ble |T] and color code symbols indicating the crystal struc¬ 
tures in Fig. [^accordingly. Elements with strong corre¬ 
lations between W and U fluctuations also have strong 
correlations between C/ipl and U. Fig. compares the 
experimentally determined structure factor S{q) of Mg 
at the triple point with that of the IPL model 

(without using any free parameters). The agreement is 
excellent. The deviation at short q vectors, related to the 
bulk modulus, is due to the mean-field g{p) treatment of 
long-ranged attractive interactions m- 

Inspired by the seminal 1972 paper by Hoover, Young 
and Grover [34] we look for a connection between den¬ 
sity scaling exponents and crystal packings: An IPL liq¬ 
uid is known |3S1 IMl HHl SIl EZl [Z3 [ZZ] to crystallize 
into a close packed (cp) fee crystal structure when the 
IPL 7 is above 2.5, while the more open body centered 
cubic (bcc) crystal is stable at lower 7 ’s (Fig. [^. The 
IPL model exhibit a polymorphic cp-bcc transition when 
7 < 2.5 also seen for many metals [10l|78l[79]. Thus, at 
the triple point elements with 7 < 2.5 are expected [M] 
to form bcc crystals while elements with 7 > 2.5 are ex¬ 


FIG. 6. Experimental structure factor of Mg at the triple 
point [ 73 ] (green squares), compared to the IPL prediction 
(red triangles). This prediction has no free parameters: The 
exponent 7 = 3.4 of the IPL potential (Eq. [^ is given by 
experimental values EU of Cv, ctp and Kt (using Eq. 0 )’ 
and a — 1.094 was chosen to match that of the solid-liquid 
coexistence shown in Fig. [T] The blue dots are the results of 
an ab initio DFT calculation. The agreement is excellent. 


pected to form cp crystals (fee or hep). Fig. shows the 
agreement of this prediction with the experimental crys¬ 
tal structure formed at the triple point (using the DFT 
7 ’s given in Table jlj). Interestingly, the predictions fails 
for several transition metals (V, Cr, Mn, Fe, Nb, Mo, Ta, 
W and Hg), most post-transition metals (Ga, In, Sn and 
Tl) and the metalloids Si and Ge. Moreover, the IPL 
model predicts that metals that form bcc structures at 
the triple point also have a low temperature cp phase - 
this is not always the case. For all of these elements we 
find significant WU correlations (see Fig. j^. Thus, we 
conclude that scale invariance can be present even when 
the elements crystal structure cannot be accurately pre¬ 
dicted by IPL-like interactions. This also illustrates the 
importance of choosing an ab initio method for investi¬ 
gating hidden scale invariance. 


F. The melting lines follow isomorphs 

The melting line follows an isomorph to a good ap¬ 
proximation [TUI nn [55] . In Table 0 this prediction is 
validated for the period-three metals by showing that 
the melting-line scaling exponent 7 m = dlnTm/dlnp 
[551 1551 1551 15T] agrees with that of the isomorph, the 
density-scaling exponent 7 . 















G. Elevated pressure behavior of Fe and P 



FIG. 7. Phase diagram of the IPL model (Eq. ^ with 
ksT — p — 1 and g{p) — 0 ) computed with the interface 
pinning method [SD]. Green triangles are the liquid sides of 
the liquid-fcc transitions, and blue diamonds are liquid-bcc 
transitions. Full lines are guides to the eye. The vertical 
dashed line indicate the bcc-fcc transition at 7 = 2.5 | 341I77 |. 
The inset show the pair energies and a representative liquid 
radial distribution at coexistence. 


Fe shows a fairly good correlation between the virial 
and the potential energy. Interestingly, the melting line 
Tm{p) reported in Refs. [HU US] coincides to a good de¬ 
gree with an extrapolation along the isomorph starting 
from the triple point, using Tm oc p'^ with 7 = 3.6 (see 
Table |!]). 

To exemplify that the behavior of elements becomes 
simpler at high pressure we compare in Fig. [^the WU 
scatter plot of Fe at the triple point (Fig. [^a)) with sim¬ 
ulation results at the pressure 310 GPa (Fig. |^b)) corre¬ 
sponding to a pressure in Earth’s liquid outer core, and 
near the melting line of pure iron [H2 [52] ■ The correla¬ 
tion coefficient increases from R = 0.95 to i? = 0.98. Fig. 
I^c) shows that WU fluctuations are uncorrelated in the 
gas-liquid coexistence regime, which is presumably as¬ 
sociated with the formation of a gas-liquid interface in 
the simulation cell and the onset of critical fluctuations. 
Likewise, the WU correlation coefficient of the non-metal 
phosphor (P) increases at high pressures from nearly un¬ 
correlated at the triple point to R = 0.30 and R = 0.52 at 
(p = 63 GPa, T = 1880 K) and {p = 210 GPa, T = 2200 
K), respectively. 

V. CONCLUDING REMARKS 
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FIG. 8 . Periodic table where green elements has a crystal 
structure at the triple point luirTSirTgi predicted by assum¬ 
ing IPL interactions (see Sec. IVE|. The IPL prediction fails 
for red elements. Zn and Cd have been colored yellow since 
their distorted hep crystals makes the IPL prediction ambigu¬ 
ous. All colored elements have hidden scale invariance seen 
in correlated WU fluctuations (see Fig. [^ and Table [^. 


TABLE 11. Gorrelation coefficients and scaling exponents for 
the period-three metals in the liquid phase at the triple point. 
Gomparing the last two columns show that the melting-line 
scaling exponent 7 ^ agrees with the density-scaling exponent 

7- 


Element 

R 

7 

7m 

Na 

0.83(0.02)“ 

1.9(0.1) 

1.7(0.5) 

Mg 

0.93(0.04) 

2.8(0.3) 

2.6(0.5) 

A1 

0.86(0.06) 

4.0(0.4) 

4.6(0.8) 


“Numbers in parenthesis indicate the statistical 
uncertainties. 


A. The Griineisen equation of state 


It has been known for a long time that at pressures high 
enough to result in non-negligible compression, solids and 
liquids generally obey the Griineisen equation of state 
from 1912 (also referred to as the Mie-Griineisen equa¬ 
tion) [HU EH- This expresses proportionality between 
pressure and energy E per volume V as follows: 


P = 1g{.p)E/V + C{p) 


(9) 


in which 7 g(p) is the Griineisen parameter (Sec. IVD) 
and C{p) the “cold pressure”, both of which are functions 
only of the density. The Griineisen equation has been ap¬ 
plied for describing condensed matter in a wide variety 
of high-pressure situations, ranging from the core of the 
Earth [52] to various forms of explosions [HU |HS]. The 
well-proven Hugoniot shock-adiabatic method is avail¬ 
able for determining 70 (p) experimentally [HU IHHl 157] . 
At high pressure in the dense fluid phase not too far 
from the melting line [HH] the virial term dominates. For 
instance, the ideal-gas contribution to the pressure is be¬ 
low five percent throughout the liquid outer core of the 
Earth. The configuration-space analog of the Griineisen 
equation is the relation 


^ — l{p)U + C{p)V. 


( 10 ) 


Equation (10) follows from the hidden scale invariance 
m- The scaling exponent j{p) has a simple relation to 
7 g(p) given in Eq. 0 m- 
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Fe at the triple point 



Potential energy [eV] 


Fe at Earth’s core pressure 
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FIG. 9. Virial potential-energy scatter plots for iron. Panel (a) shows a scatter plot of potential energy and virial of 125 iron 
atoms at the DFT triple point. Panel (b) shows that the correlation increases at elevated pressure, here chosen corresponding 
to a pressure at the Earth’s core and near the melting line of pure iron (310 GPa, 9000 K) |82l I83| . Thus iron becomes simpler 
at higher pressures. Panel (c) shows that the correlation is low in the gas-liquid coexistence region. 
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FIG. 10. High-pressure calculations of phosphorus (P) show 
that the correlation coefficient increases with increasing pres¬ 
sure. 


B. Empirical melting and freezing rules 


second-order liquid configurational entropy at crystal¬ 
lization |97j . Connecting the melting and freezing lines 
is the rule of invariant constant-volume melting entropy 

[Siin]. 
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As a consequence of our finding, many empirical melt¬ 
ing and freezing rules now find a concise explanation. 
Specihcally, a number of invariants along the melting 
line of metals (and model systems) have been known 
for years with no good explanations. These rules fol¬ 
low from hidden scale invariance, because the melting 
and freezing lines are both isomorphs m and the rules 
all involve isomorph invariants. A famous melting rule 
is the Lindemann criterion, according to which a crys¬ 
tal melts when the thermal vibrational atomic displace¬ 
ment is about 10% of the crystal’s interatomic distance 
[HHl [RnU?TT] . While our work does not imply a universal 
value of 10%, it does imply invariance of the Lindemann 
quantity along the melting line. There are also other em¬ 
pirically well-established freezing rules of invariance, for 
instance the Hansen-Verlet rule that a liquid crystallizes 
when the first peak of the structure factor reaches the 
value 2.85 [HI], the Andrade equation predicting constant 
reduced-unit viscosity along the freezing line [HSl IMj . 
the Raveche-Mountain-Streett criterion [95] of a quasiu- 
niversal ratio between maximum and minimum of the ra¬ 
dial distribution function at freezing, Lyapunov-exponent 
based criteria [96] . or the criterion of zero higher-than- 


Appendix A: Estimations of DFT triple points 

The main focus of the paper is to determine hidden 
scale invariances of elements in the low-pressure part of 
the liquid phase. An unbiased choice is to use state 
points near the DFT triple points. However, an ac¬ 
curate determination of the triple point is beyond the 
scope of the present paper, and also not relevant for 
the present work. The DFT triple point is estimated by 
performing an NpT computation at a state point close 
to, or slightly above, the experimental melting tempera¬ 
ture and ambient pressure (within an accuracy of 2 GPa 
or less, see Table |T|). The NpT ensemble is realized us¬ 
ing the Parinello-Rahman method |100j with a fictitious 
mass of 100 atomic units and a thermal coupling time 
of 0.33 ps. The finite cut-off of the plane wave basis set 
requires the exertion of an additional pressure to compen¬ 
sate for the missing basis set functions. This pressure is 
referred to as Pulay stress and it was estimated in a sepa¬ 
rate calculation. If the system froze at the experimental 
temperature, for instance due to finite size effects, the 
temperature was increased until the system stayed liquid 
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FIG. 11. Panel (a) show an example of a crystallized config¬ 
uration (Ga). The inset shows the computed diffraction pat¬ 
ten. Whenever this happened, computations was redone at 
higher temperature to avoid crystallization. Panel (b) shows 
a DFT configuration of Fe at the experimental triple point. 
DFT overestimates the triple point density by 16%; thus cav¬ 
ities are formed when iron is simulated at the experimental 
density. 


throughout the entire trajectory. An example of a crys¬ 
tallized configuration is shown on Fig. 11’a). For some 


systems, such as Ti, Ga, Sn and Hg the temperature had 
to be increased significantly compared to the experimen¬ 
tal melting temperature. However, we only found small 
variations of 7 and R upon increasing the temperature 
by a factor 2. The final state points used to compute the 
UW correlations are listed in TabU 


For Mn, Fe, and Co the computed melting density is 
25%, 16% and 15% higher than the experimental density. 
Simulations of these liquids at the experimental density 
using standard density functional theory led to internal 
surfaces and cavities between the atoms as shown in Fig. 
[ntb). The reason for this sizable deviation of the melt¬ 
ing densities from experiment needs to be investigated 
in the future. A likely explanations is strong paramag¬ 
netic fluctuations at the atomic sites that would require a 
treatment beyond standard density functional theory, for 
instance using dynamical mean field theory. Such fluc¬ 
tuations can increase the typical bond length by some 
5%. 
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